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Microscopic models of the interaction between grain boundaries (GBs) and both dislocations and 
cracks are of importance in understanding the role of microstructure in altering the mechanical 
properties of a material. A recently developed mixed atomistic and continuum method is extended to 
examine the interaction between GBs, dislocations and cracks. These calculations elucidate plausible 
microscopic mechanisms for these defect interactions and allow for the quantitative evaluation of 
critical parameters such as the stress to nucleate a dislocation at a step on a GB and the force 
needed to induce GB migration. 
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With the continuing development of more accurate, 
less expensive models for atomistic interactions and ex- 
pansion of computational resources, there is growing in- 
terest in the modeling of materials from fundamental 
principles rather than phenomenological approaches. An 
outstanding problem in this regard is the role of mi- 
crostructure in determining material properties. The in- 
fluence of microstructure (e.g. grain size and shape) on 
the mechanical properties of materials is clearly revealed, 
for example, in the yield strength and the fracture tough- 
ness [Q . A first step in the microscopic determination of 
the role of microstructure in governing such properties 
is the elucidation of plausible mechanisms whereby dis- 
locations and cracks, the primary agents of permanent 
deformation, interact with the boundaries that make up 
that microstructure. One of the key challenges posed by 
such calculations is the simultaneous operation of mul- 
tiple scales in the problem. In this letter, we present a 
recently developed model for treating multiple scales. We 
then demonstrate the application of the model to three 
examples: the deformation of a stepped grain boundary 
(GB), the interaction of lattice dislocations with GBs and 
the interaction of cracks with GBs. 

Recently, the quasicontinuum method was proposed 
to allow for a seamless treatment of multiple length 
scales . This mixed atomistic-continuum formulation 
is based on a finite element discretization of a contin- 
uum mechanics variational principle. The finite element 
method serves as the numerical engine for determining 
the energy minimizing displacement fields, while atom- 
istic analysis is used to determine the energy of a given 
configuration. This is in contrast to standard finite el- 
ement approaches where the constitutive input is made 
via phenomenological models. The method is successful 
in capturing the structure and energetics of dislocations. 

In this paper we generalize the method to allow for the 
treatment of interfaces, and show how the formulation al- 
lows for the simultaneous treatment of dislocations, ma- 
terial interfaces and cracks. We begin with the recog- 
nition that from the microscopic perspective the body 
may be regarded as a collection of N atoms, the total 
potential energy of which is given by 
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where is the position of the atom i, fi is the exter- 
nal force on that atom, and Ei is its energy as would be 
computed from an atomistic model such as the embedded 
atom method used here. One of the primary objectives 
in the formulation of the method is to eliminate the re- 
dundant atomistic degrees of freedom associated with the 
regions of the body far from extended defects and hence 
subject to displacement fields which are slowly varying 
on the atomic scale. To achieve the requisite degree of 
freedom reduction we select M representative atoms from 
the N atoms (M << N), chosen to best represent the 
energetics of the body, the positions ra{a = 1, M) of 
which serve as the reduced set of degrees of freedom. The 
body is now divided into disjoint cells such that each cell 
contains exactly one representative atom. The key ener- 
getic approximation is that the energy of all of the atoms 
in a given cell is the same as that of the cell's repre- 
sentative atom. The positions of the atoms that are not 
treated explicitly are obtained by interpolating the nodal 
values of the displacements using a finite element mesh 
which is constructed with the representative atoms as the 
nodal points. (One possible implementation of this strat- 
egy in two dimensions is to use the Voronoi polygons |^ 
surrounding the representative atoms as the cells, and 
the geometric dual of the Voronoi tiling, the Delaunay 
triangulation as the finite element mesh.) 

Given the scheme described above, the approximate 
potential energy depends only on the positions of the 
representative atoms and can be written as 
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where ria is the number of atoms represented by atom a, 
Ea is the energy of that representative atom and is 
the effective force acting on the a*'* representative atom. 
In practice, E^ is computed in two different ways. When 
the representative atom is located in a region undergo- 
ing strongly non-uniform deformation, Ea is computed 
using the usual atomistic rule in which a given atom is 
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surrounded by its complement of neighbors and the re- 
sulting energy per atom is computed. On the other hand, 
if the representative atom experiences a slowly varying 
deformation the energy is still computed atomistically, 
but with the assumption that its environment is distorted 
according to the local gradients of deformation. The de- 
tails of the criteria that dictate which scheme is used will 
be described elsewhere. The key outcome gained in the 
implementation of the strategy described above is the 
incorporation of the relevant atomistic nonlinearity and 
nonlocality that allows for the emergence of defects such 
as dislocations and cracks, without the attendant singu- 
larities that plague linear elastic analyses, or the burden 
of many redundant atomistic degrees of freedom. 

As in our earlier work with this method, a key cri- 
terion to the validity of the model is how well the re- 
sults of conventional atomistics are recovered in the con- 
text of defects of known structure. For the purposes of 
the present paper, the method must successfully repro- 
duce the known static geometric structures of GBs. As 
a test of the method, we have examined the structure of 
a range of GBs in several fee metals. For the moment, 
we have confined our attention to symmetric tilt bound- 
aries, using embedded-atom type potentials. The key 
quantitative tests of the outcome of these calculations 
are an appropriate reckoning of i) the interfacial energy, 
and ii) the interfacial structure. An indication of the 
typical energy differences between the quasicontinuum 
result and the associated direct atomistic calculation is 
demonstrated by a S5(210) GB in Au where the energy 
as obtained by conventional atomistics (676mJ/m^) |^] 
and the method described here (filOmJ/m^) are in close 
agreement. Similarly, in all of the cases we have consid- 
ered (i.e. E5(210) in Al and Cu, S3(lll) in Al, 1199(557) 
in Al, and S21(241) in Al, Au and Ni), the atomic level 
geometry at the interface obtained using the quasicontin- 
uum method advocated here was for practical purposes 
identical to that obtained using direct atomistic simula- 
tion. Consequently, it is of interest to turn the method 
to the analysis of interfacial deformation. 




FIG. 1. Step motion under applied shear stress. Atomic 
positions associated with stepped twin boundary, illustrating 
the motion of the boundary perpendicular to the boundary 
plane as a result of an applied shear stress, (a) Initial config- 
uration, (b) Final configuration after the application of the 
critical stress. 

As an example of such deformation, we consider a 
stepped twin boundary (S3(lll)) in Al that is subjected 



to a remotely applied shear stress. It is found that 
once a certain critical shear stress (w 0.033/^, where /i 
is the shear modulus relevant to the orientation of inter- 
est here) is attained, the step on the boundary serves as a 
source for ^[112] dislocations which are swept out along 
the boundary plane. The consequence of the passing of 
these dislocations is the net downward motion of the twin 
boundary shown in fig. 1. One question of importance 
concerning the type of simulations described above is that 
of the convergence of the results with respect to both sys- 
tem size and mesh geometry. We have found that altering 
the system size reveals that the critical stress of « 0.033/i 
necessary to move the boundary is reproduced to within 
±0.003/i from one simulation to the next. It is interest- 
ing to contrast this stress with a typical Peierls stress 
for a straight dislocation. For example, in the case of a 
screw dislocation in this metal, we have found that 
the Peierls stress is 0.00068/i, nearly fifty times smaller 
than the critical stress for advancing the twin boundary. 
As another comparison to set the scale of the stresses 
determined here, the stress to induce motion of the twin 
boundary can be compared with that to operate a Frank- 
Read source which is cr w fJ'b/ L, where L is the width of 
the source Q. In light of this estimate, the stress to in- 
duce motion of the twin boundary is of the same order 
as that to operate a Frank-Read source of width « 356 
(where 6 is a typical Burgers vector). Although typical 
Frank-Read sources are larger than 356 and hence oper- 
ate at even lower stresses, the stress found to stimulate 
motion of the twin boundary is still significantly smaller 
than the ideal shear strength, and is another example of 
the "lubricating" effect of heterogeneities in the motion 
of extended defects. 

Despite the existence of useful continuum models of 
dislocation-GB interactions, it remains a crucial chal- 
lenge to uncover the microscopic processes that transpire 
once the dislocation core is in the proximity of a GB. Our 
earlier work on simulating nanoindentation suggests 
the possibility of using nanoindentation induced disloca- 
tions to probe the interaction between dislocations and a 
GB. As a model system, we consider a block oriented such 
that (111) planes are positioned to allow for the emer- 
gence of dislocations which then travel to the 5]21(241) 
GB which waits approximately 200A beneath the surface 
(c.f. figure 2a). 

Because of the relatively high stacking fault energy as- 
sociated with the EAM potentials for Al we used 0, the 
dislocations nucleated at the free surface as a result of the 
indentation process are produced as rather closely spaced 
(15A) Shockley partials. As seen in the left frame in part 
(b) of fig. 2, the Shockley partials have been absorbed at 
the GB with the creation of a step at the GB. This ge- 
ometry can be rationalized on the basis of the underlying 
displacement shift complete (DSC) lattice associated 
with this symmetric GB. We find that the lattice dislo- 
cation ^[IlO] can be split into two DSC lattice vectors. 
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^[110]= ||[312] +^[241], (3) 

GBDislocation Step 

where [312] is the Burgers vector of a GB dislocation 
parallel to the GB and ^[241] is the vector associated 
with the step. 

As the load is increased, a second pair of Shockley par- 
tials is nucleated and they are not immediately absorbed 
into the GB and consequently form a pile-up (c.f. fig. 2). 
These dislocations are not absorbed until a much higher 
load level is attained. Due to the lack of dislocation ac- 
tivity in the neighboring grain, it may be concluded that 
this GB does not aid slip transmission. 
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FIG. 2. (a) Finite element mesh used to model disloca- 
tion/GB interaction. The surface marked AB is rigidly in- 
dented in order to generate dislocations at A. (b) Snapshots 
of atomic positions at different stages in the deformation his- 
tory. Absorption of the first pair of dislocations at the GB 
results in a step, while the second pair form a pile-up. 

As a final example of the synthetic view of extended 

defects afforded by this method, we consider the interac- 
tion between a brittle crack and a GB. The interaction 



of cracks and interfaces pose a variety of challenging and 
important problems. One issue that can be considered 
within the confines of the method presented here is that 
of the interaction of a crack propagating by cleavage as it 
impinges upon a GB in its path. The issues that attend 
the use of the method for considering fracture in general 
will be presented elsewhere, while here we will note the 
key elements in carrying out such simulations. 

In order to investigate the interaction between an ad- 
vancing crack and a GB, we consider the S21(241) GB 
in fee nickel. A crack is initiated in one of the grains 
by removing a single (111) plane, such that the crack tip 
is located about 2 nm from the GB. The crack is then 
loaded by applying the isotropic linear elastic displace- 
ment fields for a sharp crack at the mesh boundaries. The 
load is incrementally increased by scaling the boundary 
node displacements and allowing the interior nodes to 
relax to their minimum energy configuration. 
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FIG. 3. Crack/GB interaction. Snapshots of crack tip re- 
gion showing motion of crack tip, dislocations and GB. 



Two snapshots of the solution are shown in fig. 3. We 
show the atoms associated with that part of our finite 
element mesh that is fully refined to the atomic scale in 
the immediate vicinity of the crack tip. The surrounding 
mesh, which extends about 300 nm in each direction, has 
been removed for clarity. The dots are atomic positions, 
while the contours reveal displacement jumps across ac- 
tive slip planes, indicating the presence of dislocations. 
Frame (a) of this figure is the configuration after 4 load 
steps. The atom labeled "ct" indicates the initial loca- 
tion of the crack tip, and one can see that the crack has 
begun to propagate towards the GB by cleavage. Light 
grey slip traces emanating from the GB, such as those 
labeled "dl" and "d2" show where the stressed GB has 
emitted dislocations, some of which have been pinned 
against the edge of the fully refined region. The dashed 
line running diagonally through the figure indicates the 
initial location of the GB which moves as a result of the 
high stresses in the crack tip region. This motion is ac- 
commodated by the structural rearrangement of atoms 
in the left hand grain to lattice sites of the right hand 
grain due to shearing along atomic planes. The solid line 
through the figure indicates the location of the GB after 
migration. 

In frame (b) of figure 3, the solution after another few 
load steps is depicted. Here, the crack has reached the 
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GB and has been blunted when atoms above the plane of 
the crack again underwent a shearing deformation. This 
time, however, the right hand grain shears to match the 
structure of the left hand grain. The two straight solid 
lines indicate the new location of the GB. The result of 
this crack blunting is a significant reduction in the stress 
levels above the crack, and dislocations such as "dl" and 
"d2" in the first frame have moved back to be re-absorbed 
by the GB. Further loading of the crack leads to a contin- 
ued crack blunting due to shearing of atomic planes along 
the GB. We have studied other GBs (l2j, where the crack 
has defiected and continued to propagate along the GB, 
in contrast to the crack blunting mechanism described 
here. 
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FIG. 4. (a) Driving force (normalized by the elastic con- 
stant, cii ) as a function of position s along the GB (normal- 
ized by the lattice constant ao). (b) The same force superim- 
posed on the GB for comparison. 



To understand these results, we turn to continuum me- 
chanics which provides the basis for evaluating the en- 
ergetic origins of GB migration. Such reasoning asserts 
that the driving force on an interface is given by the jump 
in the Eshelby tensor |l3| across the interface, with this 
tensor defined as 

Pij = WSij - UkjCTkj . (4) 

W is the strain energy density, Ui^k is the k*^ component 
of the gradient in the i*^ component of displacement and 
akj is the stress tensor. Within the confines of linear 
elasticity, we have computed the driving force on the in- 
terface by using a conventional anisotropic linear elastic 
constitutive model in conjunction with the standard fi- 
nite element method to obtain the fields associated with 
the crack/GB geometry described above. Once these 
fields are obtained, the resulting driving force may be 
obtained by computing the jump in the Eshelby tensor 
across the interface. If we further assume that the GB 
migration is proportional to the driving force, the driving 
force profile may be compared directly with the bowed 
out geometry as shown in figure 4. 

In this letter, we have shown how our mixed atom- 
istic and continuum analyses has been adapted to the 
treatment of interfacial deformation. Such calculations 
demanded generalization of the original quasicontinuum 
formulation to allow for the existence of more than one 
grain at the same time. As validation of the method, we 
have computed the structure and energetics of a series of 



different GBs and found entirely satisfactory correspon- 
dence between these calculations and those resulting from 
direct atomistics. The method has been applied to three 
distinct problems: deformation of a stepped twin bound- 
ary, the interaction between dislocations and a GB, and 
the propagation of a crack into a GB. The calculations 
on the stepped GB allowed for quantitative evaluation 
of the stress to move the GB, while calculations on the 
crack/GB interaction revealed stress induced GB motion 
which can be rationalized in terms of the driving force 
on that interface as implied by the jump in the Eshelby 
tensor. The advantage of the model presented here over 
standard atomistic calculations is the significant reduc- 
tion in the computational effort through careful reduc- 
tion of the degrees of freedom. For example, the number 
of degrees of freedom associated with the mesh of fig- 
ure 2a is about 10^, while the same atomistic calculation 
would have required more than lO'' degrees of freedom. 
This approach allows for the simultaneous treatment of 
defects occurring over many length scales, ranging from 
individual dislocations to GBs and cracks. 
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